using Peridynamics

function sphere_impact(; ns=6, np=3, path)
    Ø = 0.15
    ΔX0_sphere = Ø / ns
    pos_sphere, vol_sphere = uniform_sphere(Ø, ΔX0_sphere; center_z=Ø / 2 + ΔX0_sphere)
    sphere = Body(BBMaterial{EnergySurfaceCorrection}(), pos_sphere, vol_sphere)
    failure_permit!(sphere, false)
    material!(sphere; horizon=3.1ΔX0_sphere, E=210e9, rho=8000, Gc=2000)
    velocity_ic!(sphere, :all_points, :z, -20)
    lxy, lz = 2.0, 0.1
    ΔX0_plate = lz / np
    pos_plate, vol_plate = uniform_box(lxy, lxy, lz, ΔX0_plate; center_z=-lz / 2)
    plate = Body(BBMaterial{EnergySurfaceCorrection}(), pos_plate, vol_plate)
    material!(plate; horizon=3.1ΔX0_plate, E=27e9, rho=2700, Gc=10)
    ms = MultibodySetup(:sphere => sphere, :plate => plate)
    contact!(ms, :sphere, :plate; radius=min(ΔX0_sphere, ΔX0_plate))
    vv = VelocityVerlet(steps=2000)
    job = Job(ms, vv; path=path)
    @time submit(job)
    return nothing
end

##
path = "results/benchmarks/sphere_impact"
rm(path; recursive=true, force=true)
sphere_impact(; path) # compilation
sphere_impact(; ns=20, np=8, path) # work

## Peridynamics.jl v0.3.0 - everywhere using Polyester.@batch
# compilation:
# 89.256836 seconds (47.14 M allocations: 14.914 GiB, 0.66% gc time, 7.14% compilation time)
# work:
# 45.514108 seconds (40.16 M allocations: 128.786 GiB, 8.94% gc time)

## Peridynamics.jl v0.3.1-DEV - everywhere using Threads.@threads :static
## PointNeighbors.jl v0.4.5: update_grid!(...; parallelization_backend=ThreadsStaticBackend())
# compilation:
# 106.717161 seconds (54.34 M allocations: 15.806 GiB, 0.61% gc time, 5.13% compilation time)
# work:
# 39.184317 seconds (48.59 M allocations: 129.752 GiB, 9.53% gc time)